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ABSTRACT 

The Euclidean Median (EM) of a set of points Q in an Euclidean 
space is the point x minimizing the (weighted) sum of the Euclidean 
distances of x to the points in Q. While there exits no closed-form 
expression for the EM, it can nevertheless be computed using iterative 
methods such as the Weiszfeld algorithm. The EM has classically 
been used as a robust estimator of centrality for multivariate data. It 
was recently demonstrated that the EM can be used to perform robust 
patch-based denoising of images by generalizing the popular Non- 
Local Means algorithm. In this paper, we propose a novel algorithm 
for computing the EM (and its box-constrained counterpart) using 
variable splitting and the method of augmented Lagrangian. The 
attractive feature of this approach is that the subproblems involved in 
the ADMM-based optimization of the augmented Lagrangian can be 
resolved using simple closed-form projections. The proposed ADMM 
solver is used for robust patch-based image denoising and is shown 
to exhibit faster convergence compared to an existing solver. 

Index Terms — Image denoising, patch-based algorithm, robust¬ 
ness, Euclidean median, variable splitting, augmented Lagrangian, 
alternating direction method of multipliers (ADMM), convergence. 

1. INTRODUCTION 

Suppose we are given a set of points ai,..., an e R'*, and some 
non-negative weights wi^... ^Wnio these points. Consider the prob¬ 
lem 

n 

min Wk \\x - ak \\2 (1) 

£cGC ^^ 
k=l 

where || • ||2 is the Euclidean norm and C C is closed and convex. 
This is a convex optimization problem. The (global) minimizer of 

(1) when C is the entire space is called the Euclidean median (also 
referred to as the geometric median) [1]. The Euclidean median has 
classically been used as a robust estimator of centrality for multivari¬ 
ate data [1]. This multivariate generalization of the scalar median also 
comes up in transport engineering, where (1) is used to model the 
problem of locating a facility that minimizes the cost of transportation 

[2] . More recently, it was demonstrated in [3, 4] that the Euclidean 
median and its non-convex extensions can be used to perform robust 
patch-based regression for image denoising by generalizing the poplar 
Non-Local Means algorithm [5, 6]. 

Note that the minimizer of (1) is simply the projection of the 
Euclidean median onto C. Unfortunately, there is no simple closed- 
form expression for the Euclidean median, even when it is unique [7]. 
Nevertheless, the Euclidean median can be computed using numerical 
methods such as the iteratively reweighted least squares (IRLS). One 
such iterative method is the so-called Weiszfeld algorithm [8]. This 
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particular algorithm, however, is known to be prone to convergence 
problems [2]. Geometric optimization-based algorithms have also be 
proposed for approximating the Euclidean median [7]. More recently, 
motivated by the use of IRLS for sparse signal processing and Li 
minimization [9, 10], an IRLS-based algorithm for approximating 
the Euclidean median was proposed in [3]. More precisely, here the 
authors consider a smooth convex surrogate of (1), namely, 

n 

(11® - Ofelli(e > 0), (2) 

k = l 

and describe an iterative method for optimizing (2). The convergence 
properties of this iterative algorithm was later studied in [11]. 

In the last few years, there has been a lot of renewed interest 
in the use of variable splitting, the method of multipliers, and the 
alternating direction method of multipliers (ADMM) [12] for various 
non-smooth optimization problems in signal processing [13, 14, 15] . 
Motivated by this line of work, we introduce a novel algorithm for 
approximating the solution of (1) using variable splitting and the aug¬ 
mented Lagrangian in Section 2. An important technical distinction 
of the present approach is that, instead of working with the smooth 
surrogate (2), we directly address the original non-smooth objective 
in (1). The attractive feature here is that the subproblems involved 
in the ADMM-based optimization of the augmented Lagrangian can 
be resolved using simple closed-form projections. In Section 3, the 
proposed algorithm is used for robust patch-based denoising of im¬ 
ages following the proposal in [3]. In particular, we incorporate the 
information that the pixels of the original (clean) patches take val¬ 
ues in a given dynamic range (for example, in [0,255] for grayscale 
images) using an appropriately defined C. Numerical experiments 
show that the iterates of the proposed algorithm converge much more 
rapidly than the IRLS iterates for the denoising method in question. 
Moreover, the reconstruction quality of the denoised image obtained 
after few ADMM iterations is often substantially better than the IRLS 
counterpart. 

2. EUCLIDEAN MEDIAN USING ADMM 

Since we will primarily work with (1) in this paper, we will overload 
terminology and refer to the minimizer of (1) as the Euclidean me¬ 
dian. The main idea behind the ADMM algorithm is to decompose 
the objective in (1) into a sum of independent objectives, which can 
then be optimized separately. This is precisely achieved using vari¬ 
able splitting [13, 14]. In particular, note that we can equivalently 
formulate (1) as 

n 

min Lc(z) + fk{xk) 

( 3 ) 

subject to Xk — z = 0 (A; = 1,..., n), 




where 


where 


(4) 

and ic is the indicator function of C [13]. In other words, we artifi¬ 
cially introduce one local variable for each term of the objective, and 
a common global variable that forces the local variables to be equal 
[16]. The advantage of this decomposition will be evident shortly. 

Note that in the process of decomposing the objective, we have 
introduced additional constraints. While one could potentially use 
any of the existing constrained optimization methods to address (3), 
we will adopt the method of augmented Lagrangian which (among 
other things) is particularly tailored to the separable structure of the 
objective. The augmented Lagrangian for (3) can be written as 

n 

£ = ic{z) + ^fk{xk) + vlixk - z) + ^\\xk - «||i| (5) 

k = l 

where ... are the Lagrange multipliers corresponding to the 
equality constraints in (3), and /x > 0 is a penalty parameter [12, 13]. 

We next use the machinery of generalized ADMM to optimize C 
jointly with respect to the variables and the Lagrange multipliers. The 
generalized ADMM algorithm proceeds as follows [13, 16]: In each 
iteration, we sequentially minimize C with respect to the variables 
cci, CC 2 , • • •, ccn, ^ in a Gauss-Seidel fashion (keeping the remain¬ 
ing variables fixed), and then we update the Lagrange multipliers 
y^^... ^y^ using a fixed dual-ascent rule. In particular, note that the 
minimization over cc^ at the (t + l)-th iteration is given by 

= argmin fk{x) + {x - ^ II® - II 2 (6) 

X ^ 


/(*) = A 11® - u||2. (10) 

Note that the objective in (9) is closed, proper, and strongly convex, 
and hence cc* exists and is unique. In fact, we immediately recognize 
cc* to be the proximal map of / evaluated at it [13, 17]. 

Definition 2.1 For any f : (—(X), od] that is closed, proper, 

and convex, the proximal map R^ is defined to be 

= argmin/(®) + f I® - v||i {vGlC). (11) 

X ^ 

As per this definition, cc* = ^/(v). 

Note that (10) involves the L 2 norm, which unlike the Li norm 
is not separable. A coordinate-wise minimization of (9) is thus not 
possible (as is the case for the Li norm leading to the well-known 
shrinkage function). However, we have the following result on proxi¬ 
mal maps [13, 17]. 

Theorem 2.2 Let f : R^ (— (X), (X)] be closed, proper, and con¬ 
vex. Then we can decompose any x G R^ as 

cc = T^/(cc) + (cc), (12) 

where /* : R^ 1 -^ (— (X), (X)] is the convex conjugate of f, 

f*{y) = max x'^y - f{x) (y e R^*). (13) 

X 

The usefulness of (12) is that both /* and and its proximal map can 
be obtained in closed-form. Indeed, by change-of-variables. 


where denotes the variable 2 ; at the end of the t-th iteration. A 
similar notation is used to denote the Lagrange multipliers at the end 
of the t-th iteration. 

Next, by combining the quadratic and linear terms and discarding 
terms that do not depend on 2 ;, we can express the partial minimiza¬ 
tion over 2 ; as 

N 

= argmin lc{z) + | ^ Ik - 

^ ^ k = l 

= argmin ||^-4*+')||i 
zee 

where 

^ k=l 

Therefore, 

where (^c) denotes the projection of x onto the convex set C. Note 
that we have used the most recent Xk variables in (7). Finally, we 
update the Lagrange multipliers using the following standard rule: 

(8) 

We loop over (6), (7) and (8) until some convergence criteria is 
met, or up to a fixed number of iterations. For a detailed exposition 
on the augmented Lagrangian, the method of multipliers, and the 
ADMM, we refer the interested reader to [12, 13, 16]. 

We now focus on the update in (6). By combining the linear and 
quadratic terms (and discarding terms that do not depend on x), we 
can express (6) in the following general form: 


/*(y) = max x^y - A |1® - u ||2 

X 

= y + max x^y — A ||a3||2. 

X 

Now, if ||y ||2 < A, then by Cauchy-Schwarz, x^y < A ||cc|| 2 . In 
this case, the maximum over x is 0. On the other hand, if ||i /||2 > A, 
then setting cc = (t > 0), we have 

x^y - A ||a3||2 = t \\y \\2 i\\y \\2 - A) > 0, 


which can be made arbitrarily large by letting t ^ 00 . Thus, 

f*(y) = u^y + tisiy), (14) 


where lb is the indicator function of the ball B with centre 0 and 
radius A, 


I'siy) = 



if l|y||2 < A, 
else. 


Having obtained (14), we have from (11), 


(v) = argmin x lb{x)-\\ x 

X 2 

= argmin L® - {v - u)\\l, 
xeB 2 


v|| 


2 

2 


which is precisely the projection of v — it onto B. This is explicitly 
given by 


= min(A, ||v-m|| 2 ) (15) 

Combining (12) with (15), we have 


* 


argmin /(cc) + 

X 



(9) ’I'/k) = 'I'a(v,'w) = v - min(A, ||v - «|| 2 ) -ji - rp 

Ik - m||2 


X 


(16) 






It is reassuring to note that tt) equals v when A = 0, and 

equals u when A = oo. In the context of the ADMM update (6), note 
that 

X =u = ttk, mid V = 

The overall ADMM algorithm (called EM-ADMM) for computing the 
Euclidean median is summarized in Algorithm 1. 


Data: Dimension d, points ai,..., an G R^, and fi>0. 

Result: z. 

Initialize: yj,..., € R'^; 

loop 

r = 0; 

for /c = 1,..., n do 

r = r + Xk + 

end 

for /c = 1,..., n do 

I j/fe = y/c + 

end 

end 

Algorithm 1: Euclidean Median using ADMM (EM-ADMM) 

We next demonstrate how EM-ADMM can be used for robust 
patch-based denoising of images. We also investigate its convergence 
behavior in relation to the IRES algorithm in [3]. 

3. APPLICATION: ROBUST PATCH REGRESSION 

In the last ten years or so, some very effective patch-based algorithms 
for image denoising have been proposed [5, 18, 19, 20, 21]. One of 
the outstanding proposals in this area is the Non-Local Means (NLM) 
algorithm [5]. Several improvements and adaptations of NLM have 
been proposed since its inception, some of which currently offer state- 
of-the-art denoising results [21]. Recently, it was demonstrated in [3] 
that the denoising performance of NLM can be further improved by 
introducing the robust Euclidean median into the NLM framework. 

We consider the denoising setup where we are given a noisy 
image g — {gi)iei {I is some linear ordering of the pixels in the 
image) that is derived from some clean image / = (/via 

= + 

where ^ = {^i)iei is iid A/'(0,1). The denoising problem is one of 
estimating the unknown f from the noisy measurement g. 

Eor any pixel z G /, let patch denote the restriction of gf to a 
square window centered at i. Letting k be the length of this window, 
this associates every pixel i with a patch-vector of length k^. For 
z, j G /, define the weights Wij = exp(—||Pi — Pjlj^/Zz^), where 
/z > 0 is a smoothing parameter. For every z G /, we consider the 
following optimization problem 

Pi = argmin tc(P) + Wij ||P - Pj || 2 . (17) 

jes(i) 

Here S{i) denotes the neighborhood pixels of i, namely, those pixels 
that are contained in a window of size S x S centred at z. The convex 
set C is defined to be the collection of patches whose pixel intensities 
are in the dynamic range [/,zz]; e.g., I = 0 and u = 255 for a 



Fig. 1. PSNR evolution with iterations for the Barbara image at two 
different noise levels. For the ADMM, we used fi — 10and for 
the IRLS, we set £ = 10In either case, we use the noisy patch to 
seed the iterations (please see main text for further details). 


grayscale image. The projection onto C can be computed coordinate- 
wise at negligible cost. In particular, for 1 < z < d, 

{ x[i\ if I < x[i\ < u, 

I if x[i] < I, and 

u if x[i] > u. 

Notice that (17) is exactly the optimization in (1). In other words, 
we denoise each patch by taking the Euclidean median of its neigh¬ 
bouring noisy patches. The center pixel of the denoised patch is taken 
to be the denoised pixel. The overall denoising scheme is called the 
Non-Local Euclidean Medians (in short, NLEM) [3]. 


Table 1. Denoising performance of (a) NLM, (b) NLEM-IRLS and 
(c) NLEM-ADMM at cr = 10, 20,40,60, and 80 (PSNRs averaged 
over 10 noise realizations). Parameters: iS = 21, /c = 7, /z = lOcr. 


Image 

Method 



PSNR (dB) 




(a) 

34.22 

29.78 

25.20 

23.37 

22.35 

House 

(b) 

32.66 

29.01 

26.68 

24.78 

23.46 


(c) 

34.13 

30.77 

27.04 

24.87 

23.47 


(a) 

33.24 

29.31 

26.17 

24.54 

23.64 

Lena 

(b) 

32.54 

29.48 

27.31 

25.38 

24.37 


(c) 

33.57 

30.38 

27.59 

25.42 

24.38 


(a) 

32.32 

27.66 

23.11 

21.03 

19.99 

Peppers 

(b) 

30.95 

26.95 

24.31 

22.84 

21.24 


(c) 

32.14 

28.54 

25.06 

22.98 

21.26 


(a) 

32.37 

27.39 

23.53 

22.05 

21.34 

Barbara 

(b) 

30.77 

27.28 

25.55 

23.77 

22.61 


(c) 

32.43 

28.84 

25.67 

23.77 

22.62 


The denoising experiments were performed on several grayscale 
test images. We used the standard NLM parameters for all the experi¬ 
ments [5]: S — 21, /c = 7, and h = lOcr. The proposed EM-ADMM 
algorithm is used for computing the solution of (17). Notice that the 




























main computation here is that of determining the distances in (16). 
Since we require a total of such distance evaluations in k‘^ dimen¬ 
sions, the complexity is 0(5'^/c^) per iteration per pixel. Incidentally, 
this is also the complexity of the IRLS algorithm in [3]. For all the 
denoising experiments, we initialized the Lagrange multipliers in 
EM-AD MM to zero vectors. We considered two possible initializations 
for namely, the noisy patch and the patch obtained using NLM. 
Exhaustive experiments showed that the best convergence results are 
obtained for /x ~ 10“^. For example, figure 1 shows the evolution 
of the peak-signal-to-noise ratio (PSNR) over the first ten iterations 
obtained using the ADMM algorithm. Also show in this figure is the 
evolution of the PSNR for the IRLS algorithm from [3]. In either 
case, we used the noisy patch as the initialization. Notice that the 
increase in PSNR is much more rapid with ADMM compared to 
IRLS. In fact, in about 4 iterations, the optimal PSNR is attained with 
ADMM. Notice that at cr = 10, the PSNR from IRLS grows much 
more slowly compared to ADMM. 



3.2 


123456789 10 

iterations 

Fig. 2. Convergence of EM-ADMM and IRLS at a given pixel of the 
House image at cr = 40. The inset shows the position of the particular 
pixel where the patch regression is performed. The objective functions 
shown in the plot correspond to (1) and (2) that are respectively 
optimized by ADMM and IRLS (note that the objectives are very 
close since we use £ = 10“® in (2)). The noisy image patch was used 
to initialize both algorithms. For EM-ADMM, the objective function 
converges to the optimal value (up to to three decimal places) in just 
2 iterations. The convergence is relatively slow for IRLS (accuracy 
up to three decimals obtained after 6 iterations). 

The above observation is also supported by the convergence result 
shown in figure 2 for the House image. Notice that ADMM converges 
in just 2 iterations, while IRLS takes about 6 iterations to attain the 
same accuracy (we recall that the cost per iteration is comparable). 
In general, we observed that the within 2 iterations the ADMM result 
converges to an accuracy that is sufficient for the present denoising 
application. At higher noise levels (cr > 60), about 4 iterations are 
required. This applies to all the test images that we have experimented 
with. 

The denoising results for some test images obtained using NLM 
and NLEM are provided in Table 1. For NLEM, we used both the 
ADMM and IRLS solvers. Following the previous observations, we 
used 4 iterations for both solvers. For a < 60, we used the noisy 


patch to initialize the iterations, and for a > 60 we used the NLM 
patch. This initialization scheme was found to give the best PSNR 
results. At large noise levels, there is not much difference after 4 
iterations. However, at low noise levels, the PSNR obtained using 
ADMM is often substantially larger than that obtained using IRLS 
after 4 iterations. The reason for this is that IRLS converges really 
slow in such situations. This is evident from the PSNR plots in figure 
1 at cr = 10. For a visual comparison, a particular denoising result 
for the Barbara image obtained using NLM and NLEM (using the 
ADMM solver) is shown in figure 3. Roughly speaking, at the cost 
of just four NLMs, we are able to improve the PSNR by more than 2 
dB. Notice that the NLEM result is much more sharp compared to 
the NLM counterpart. 


(a) Barbara (512 x 512). (b) Corrupted (PSNR =16.15 dB) 


(c) NLM (PSNR =23.53 dB) (d) NELM (PSNR = 25.67 dB). 

Fig. 3. Denoising results using NLM and NLEM at cr = 40. We 
used the proposed EM-ADMM algorithm for computing the Euclidean 
medians at each pixel (please see the main text for other details about 
the experiment). The EM-ADMM algorithm was initialized using the 
noisy patch and run for four iterations. 


4. CONCLUSION 

We proposed a new algorithm for computing the (constrained) Eu¬ 
clidean median using variable splitting and the augmented Lagrangian. 
In particular, we demonstrated how the ADMM-based optimization of 
the augmented Lagrangian can be resolved using simple closed-form 
projections. The proposed algorithm was used for image denoising 
using the Non-Local Euclidean Medians, and was generally found 
to exhibit faster convergence compared to the IRLS algorithm. One 
interesting direction that we did not pursue is to adapt the penalty fi at 
each iteration, which is known to speed up the convergence [16]. Yet 
another interesting direction is the use of accelerated ADMM [22] to 
further speed up the convergence. 
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